data <- list()

files <- dir("./data", recursive=T, full.names=T)

files.psms <- files[grep(pattern="_psms", files)]
files.proteins <- files[grep(pattern="_proteins", files)]
files.peptides <- files[grep(pattern="_peptides", files)]

data$psms <- c()
# file <- files.psms[1]
for(file in files.psms){
          if((file.info(file)[["size"]]) == 0){
                    print(paste("file", file, " is empty!!"))
                    next
          }
          
          id <- substr(file, regexpr("/", file) + 6, regexpr(pattern="_", file) -1)
          cat(sprintf("reading psms of project %s\n", id))
          file.data <- read.delim(file)
          Projectid <- rep(id, times=nrow(file.data))
          file.data <- data.frame(Projectid, file.data)
#           head(file.data)
          data$psms <- rbind(data$psms, file.data)
}

data$peptides <- c()
for(file in files.peptides){
          if((file.info(file)[["size"]]) == 0){
                    print(paste("file", file, " is empty!!"))
                    next
          }
          
          id <- substr(file, regexpr("/", file) + 6, regexpr(pattern="_", file) -1)
          cat(sprintf("reading peptides of project %s\n", id))
          file.data <- read.delim(file)
          Projectid <- rep(id, times=nrow(file.data))
          file.data <- data.frame(Projectid, file.data)
          #           head(file.data)
          data$peptides <- rbind(data$peptides, file.data)
}


data$proteins <- c()

for(file in files.proteins){
          if((file.info(file)[["size"]]) == 0){
                    print(paste("file", file, " is empty!!"))
                    next
          }
          
#           file <- files.proteins[6]
          id <- substr(file, regexpr("/", file) + 6, regexpr(pattern="_", file) -1)
          cat(sprintf("reading proteins of project %s\n", id))
          
          
          
          file.data <- read.delim(file, row.names=NULL)
          Projectid <- rep(id, times=nrow(file.data))
          file.data <- data.frame(Projectid, file.data)
          
          if(NROW(file.data) == 0){
                    next          
          }
          
          colnames(file.data) <- c("Projectid", "Protein","Equivalent.proteins"
                                   ,"Group.class","n.peptides","n.spectra","n.peptides.validated","n.spectra.validated"
                                   ,"MW","NSAF","NA","p.score","p","Decoy","Validated","Description")
          
          file.data <- head(file.data[,-c(3,11)])
          data$proteins <- rbind(data$proteins, file.data)
}


# Only keep the validated entries.
cat("filtering only validated entries")
data$psms <- data$psms[data$psms$Validated==T, ]
data$peptides <- data$peptides[data$peptides$Validated==T, ]
data$proteins <- data$proteins[data$proteins$Validated==1, ]



# Add a PSM identifier
PSM.Id <- 1:nrow(data$psms)
data$psms <- data.frame(PSM.Id, data$psms)

psm.varmods <- c()
psm.varmods.entry <- c()
for(modseq in data$psms$Variable.Modification.s.){
          
#           modseq <- "acetylation of protein n-term(1), carbamidomethyl c(16), deamidation of n and q(11), oxidation of m(15)"
          mods <- unlist(strsplit(modseq, split=","))
          curr.mods <- c()
          for(mod in mods){
                    mod <- unlist(strsplit(mod, "[(]"))
                    if(length(mod) > 1){
                              mod <- mod[1]
                              mod <- sub("^ +", "", mod)
                              curr.mods <- c(curr.mods, mod)
                              psm.varmods <- c(psm.varmods, mod)
                    }
          }
          
          if(length(curr.mods) >0){
                    psm.varmods.entry <- c(psm.varmods.entry, paste(curr.mods, collapse=";"))
          }else{
                    psm.varmods.entry <- c(psm.varmods.entry, NA)
          }
}
psm.varmods <- unique(psm.varmods)
head(psm.varmods.entry)

# Ok, now build a PTM matrix to align with the PSM rows
d <- c()
row.none <-  rep(FALSE, length(psm.varmods))

for(mods in psm.varmods.entry){
#           mods <- psm.varmods.entry[787]
          row <- row.none
          if(is.na(mods) == FALSE){
                    mods <- unlist(strsplit(mods, ";"))
                    for(mod in mods){
                              index <- grep(pattern=mod, x=psm.varmods)
                              row[index] <- TRUE
                    }
          }
          d <- rbind(d, row)
}
d <- data.frame(d, row.names=NULL)
colnames(d) <- paste("PTM", gsub(pattern=" ", replacement="_", psm.varmods), sep="_")
psm.ptms <- d

data$psms <- data.frame(data$psms, psm.ptms)
tail(data$psms)

# Extract Protein Inference state for peptides
psm.prot.inf <- c()

for(prot in data$psms$Protein.s.){
          x <- sub(" +$", "", prot)
          l <- length(unlist(strsplit(gsub(x=x, pattern=" ", replacement=","), ",")))
          psm.prot.inf <- c(psm.prot.inf, l)
}

data$psms <- data.frame(data$psms, psm.prot.inf)



## Load PRIDE experiment annotations
# Data
file <- "/Users/kennyhelsens/Proteomics/Projects/1206/1206_relims_pride/projects_120620.csv"
#get column count of each row
st.fields <- count.fields(file, sep=";")
#get max. column count
max.fields <- max(st.fields)
#define column classes
colClasses <- c(rep("character", 8),  rep("numeric", 3))
#define column names
colNames <- c("experimentId", "description", "project", "species", "taxid", "tissue", "BTO", "mods",
              "spectra", "proteins", "peptides" )

#read file
pride <- read.delim(file, sep=",", as.is=TRUE, header=FALSE, fill=TRUE, colClasses=colClasses, col.names=colNames)
# summary(pride[,c(1,5,6,8,9,10,11)])

# Make project counts
projectids <- unique(data$psms$Projectid)

psm.per.project <- data.frame(table(data$psms$Projectid))
colnames(psm.per.project) <- c("experimentId", "relims.psm.count")

pep.per.project <- data.frame(table(data$peptides$Projectid))
colnames(pep.per.project) <- c("experimentId", "relims.pep.count")

prot.per.project <- data.frame(table(data$proteins$Projectid))
colnames(prot.per.project) <- c("experimentId", "relims.prot.count")

data$counts <- merge(merge(psm.per.project, pep.per.project), prot.per.project)

data$counts <- merge(data$counts, y=pride[,c("experimentId", "spectra", "peptides", "proteins", "project", "tissue")], by="experimentId")
colnames(data$counts) <- c("experimentId", "psm.count.relims", "peptide.count.relims", "protein.count.relims",
                           "spectra.count.pride", "psm.count.pride", "protein.count.pride", "project.pride", "tissue.pride")

# Set " N\A" to NA
data$counts$tissue.pride[which(data$counts$tissue.pride == " N/A")] <- NA

min.decoy.pep <- min(data$psms[which(data$psms$Decoy==1),c("p")])
data$psms$p[which(data$psms$p <= 0.00001)] <- min.decoy.pep

save.image("./data/report_ws.RData")

